% Create Cross wavelet spectra for CBI station
% during Tohuku (2011-03-11) EQ
% for dates 2011-03-01 to 2011-03-20


%% read CBI data

fname = '/Users/manojnair/Downloads/cbi_march_2011.txt';
CBI = read_iaga2002_1min(fname);

fname = '/Users/manojnair/Downloads/mmb_2011_march.txt';
MMB = read_iaga2002_1min(fname);


%% Analize. CBI = local and MMB = remote 


local_fday = CBI.MINUTELY(:,1)';
local_h =  sqrt(CBI.MINUTELY(:,2).^2+ CBI.MINUTELY(:,3).^2)';    
local_z =  CBI.MINUTELY(:,4)'; 

remote_fday = MMB.MINUTELY(:,1)';
remote_h =  sqrt(MMB.MINUTELY(:,2).^2+ MMB.MINUTELY(:,3).^2)';    
remote_z =  MMB.MINUTELY(:,4)';




%b1 = min(local_fday):(1/24):max(local_fday);

b1 = min(local_fday):(2/24):max(local_fday);
sp=spline(b1,local_h/spline(b1,eye(length(b1)),local_fday));
local_h_f=local_h-ppval(local_fday,sp);
%clear local_h;

sp=spline(b1,local_z/spline(b1,eye(length(b1)),local_fday));
local_z_f=local_z-ppval(local_fday,sp);
%clear local_z;


sp=spline(b1,remote_h/spline(b1,eye(length(b1)),remote_fday));
remote_h_f=remote_h-ppval(remote_fday,sp);
%clear remote_h;

sp=spline(b1,remote_z/spline(b1,eye(length(b1)),local_fday));
remote_z_f=remote_z-ppval(local_fday,sp);
%clear remote_z;



%% wavelet analysis 

Sc = 1:0.5:40;

%waven = 'morl';
waven = 'cgau4';

perd = 1./scal2frq(Sc,waven,60); % 60 seconds DELTA intrvl

remotezw = cwt(remote_z_f,Sc,waven);%Complex wavelet transfor
%clear remote_z_f;

localzw = cwt(local_z_f,Sc,waven);%Complex wavelet transform
%clear local_z_f;

localhw = cwt(local_h_f,Sc,waven);%Complex wavelet transform
%clear local_h_f;

remotehw = cwt(remote_h_f,Sc,waven);%Complex wavelet transformf
%clear remote_h_f;



Hxy = localhw .* conj(remotehw);
L = abs(Hxy) > 100;
Hxy(L) = 100;
ss = sqrt(abs(Hxy));


cc = max(max(ss)) - abs(ss);
weight = (cc./max(max(cc)));

%clear cc ss L localhw remotehw ;
subplot(211);
imagesc(local_fday,perd/60,(abs(localzw)));
caxis([0,4]);
subplot(212);
imagesc(local_fday,perd/60,(abs(localzw).*weight));
caxis([0,4]);
%% plot

f = figure;
set(f,'Position',[ 34         562        1876         489]);
b(1) = min(local_fday);
b(2) = max(local_fday);

% b(1) = datenum(2010,2,27,0,0,0);
% b(2) = datenum(2010,2,28,0,0,0);

ax(1)=subplot(311);
imagesc(local_fday,perd/60,abs(localzw));
set(gca,'FontSize',16);
caxis([0,3]);

% datetick;
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Wavelet amplitude of local Z','FontSize', 18,'color','white');
set(gca,'xtick',[]);
%colorbar('peer',ax(1),[0.9083 0.7116 0.02452 0.2106],'FontSize',16);



ax(2)=subplot(312)
imagesc(local_fday,perd/60,abs(Hxy))
caxis([0,100]);
set(gca,'FontSize',16);
% b=axis;
% datetick;
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Cross wavelet amplitude between local H and remote H','FontSize', 18,'color','white');
set(gca,'xtick',[]);
%colorbar('peer',ax(2),[0.9083 0.411 0.02452 0.2106],'FontSize',16);


ax(3)=subplot(313)
imagesc(local_fday,perd/60,(abs(localzw).*weight))
set(gca,'FontSize',16);
caxis([0,3]);
datetick('x','keeplimits');
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Weighted wavelet amplitude of local Z ','FontSize', 18,'color','white');
%colorbar('peer',ax(3),[0.9083 0.1304 0.02452 0.2106],'FontSize',16);
xlabel('Date in 2011 mm/dd');
ylabel('Period in minutes');